################################################################################
## TITLE:		  figures_2spe_model_1980_2024
## AUTHORS:   P. Mongrain, R. Nadeau, B. Jérôme and V. Jérôme
## PAPER:     State-Level Forecasts for the 2024 U.S. Presidential Election
## DATE:		  December 2024


################################################################################
## INSTALL AND LOAD PACKAGES (the packages below first need to be installed)

# install.packages("dplyr")
# install.packages("effects")
# install.packages("estimatr")
# install.packages("ggeffects")
# install.packages("ggplot2")
# install.packages("ggpubr")
# install.packages("ggridges")
# install.packages("grid")
# install.packages("gridExtra")
# install.packages("haven")
# install.packages("lemon")
# install.packages("miceadds")
# install.packages("png")
# install.packages("readxl")
# install.packages("statebins")
# install.packages("statebins")
# install.packages("sjmisc")
# install.packages("tidyverse")
# install.packages("usmap")

library(dplyr)
library(effects)
library(estimatr)
library(ggeffects)
library(ggplot2)
library(ggpubr)
library(ggridges)
library(grid)
library(gridExtra)
library(haven)
library(lemon)
library(miceadds)
library(png)
library(readxl)
library(statebins)
library(sjPlot)
library(sjmisc)
library(tidyverse)
library(usmap)

# Set working directory

getwd()

setwd("") #enter file path here

################################################################################
## LOAD DATA
################################################################################

# Load dataset

df <- read_dta("data_2spe_model_1980_2024.dta")

# Define variables as factor or numeric

df$twoincv <- as.numeric(df$twoincv)
df$ptwoincv <- as.numeric(df$ptwoincv)
df$midterms2 <- as.numeric(df$midterms2)
df$legcont <- as.numeric(df$legcont)
df$unemp <- as.numeric(df$unemp)
df$jpa <- as.numeric(df$jpa)
df$incpres <- as.factor(df$incpres)
df$ppi5220 <- as.numeric(df$ppi5220)
df$ppi8020 <- as.numeric(df$ppi8020)
df$chavp <- as.numeric(df$chavp)
df$anderson <- as.factor(df$anderson)
df$perot92 <- as.factor(df$perot92)
df$perot96 <- as.factor(df$perot96)
df$dhsc <- as.factor(df$dhsc)
df$rhsc <- as.factor(df$rhsc)
df$dcds <- as.factor(df$dcds)
df$dcrs <- as.factor(df$dcrs)
df$fx <- as.factor(df$fx)

################################################################################
## EXTENDED MODEL
################################################################################

# Fit model

fit.extended <- lm_robust(twoincv ~ ptwoincv + midterms2 + legcont + unemp + jpa*incpres + ppi5220 +
                  ppi8020 + chavp + anderson + perot92 + perot96 + dhsc + rhsc + 
                  dcds + dcrs + fx, data = df, clusters = fx)

# Predicted probabilities for Job approval x Incumbent status interaction ()

inter.e <- ggpredict(model=fit.extended, terms=c("jpa", "incpres"))

inter.e.df <- as.data.frame(inter.e) #save as dataframe

inter.e.df # data used for Figure 2, panel (a), in the article

# Create interaction graph (Figure 2, panel (a), in the article)

tiff("inter_e_g.tiff", units="in", width=6, height=4, res=1200)

int.e.g <- ggplot(inter.e.df, aes(x=x, y=predicted, color=group, group=group)) +
  geom_line(size=.5, show.legend = TRUE, aes(linetype=group, color=group)) +
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high, fill=group), alpha=0.3, linetype=0, show.legend = FALSE) +
  ggtitle("(a) Extended Model") +
  scale_x_continuous(name = "Job Approval Rating", limits=c(0, 100), 
                     breaks = seq(0, 100, 20)) +
  scale_y_continuous(name = "Two-Party Vote Share", limits=c(30, 80), 
                     breaks = seq(30, 80, 10)) + theme_bw() +
  theme(axis.text.x = element_text(size = 11), axis.text.y = element_text(size = 11),
        axis.title = element_text(size = 11),
        axis.title.y = element_text(margin = margin(t = 0, r = 5, b = 0, l = 0)),
        axis.title.x = element_text(margin = margin(t = 5, r = 0, b = 0, l = 0)),
        legend.title = element_text(size = 11), legend.text=element_text(size=11),
        legend.key.size = unit(.6, 'cm'),
        #legend.position = c(0.8, 0.2),
        legend.spacing.y = unit(.01, 'lines'),
        plot.title = element_text(size=11, hjust = 0.5, colour="black"),
        #panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        plot.margin=unit(c(0,.01,.3,0), "cm")) + 
  scale_colour_manual(values=c("#F8766D", "#00BFC4"), name="Incumbent Running?", 
                      labels=c("No","Yes"), guide = guide_legend(title.position = "top",
                                                                                   reverse=TRUE)) +
  scale_fill_manual(values=c("#F8766D", "#00BFC4"), name="Incumbent Running?", 
                    labels=c("No","Yes"), guide = guide_legend(title.position = "top",
                                                                                 reverse=TRUE)) +
  scale_linetype_manual(values=c("longdash","solid"), name="Incumbent Running?",
                        labels=c("No","Yes"), guide = guide_legend(title.position = "top",
                                                                                     reverse=TRUE))

int.e.g

dev.off()

################################################################################
## SIMPLIFIED MODEL
################################################################################

# Fit model

fit.simplified <- lm_robust(twoincv ~ ptwoincv + jpa*incpres + ppi5220 +
                            ppi8020 + chavp + fx, data = df, clusters = fx)

# Predicted probabilities for Job approval x Incumbent status interaction

inter.s <- ggpredict(model=fit.simplified, terms=c("jpa", "incpres"))

inter.s.df <- as.data.frame(inter.s) #save as dataframe

inter.s.df # data used for Figure 2, panel (b), in the article

# Create interaction graph (Figure 2, panel (b), in the article)

tiff("inter_s_g.tiff", units="in", width=6, height=4, res=1200)

int.s.g <- ggplot(inter.s.df, aes(x=x, y=predicted, color=group, group=group)) +
  geom_line(size=.5, show.legend = TRUE, aes(linetype=group, color=group)) +
  geom_ribbon(aes(ymin=conf.low, ymax=conf.high, fill=group), alpha=0.3, linetype=0, show.legend = FALSE) +
  ggtitle("(b) Simplified Model") +
  scale_x_continuous(name = "Job Approval Rating", limits=c(0, 100), 
                     breaks = seq(0, 100, 20)) +
  scale_y_continuous(name = "Two-Party Vote Share", limits=c(30, 80), 
                     breaks = seq(30, 80, 10)) + theme_bw() +
  theme(axis.text.x = element_text(size = 11), axis.text.y = element_text(size = 11),
        axis.title = element_text(size = 11),
        axis.title.y = element_text(color="white", margin = margin(t = 0, r = 5, b = 0, l = 0)),
        axis.title.x = element_text(margin = margin(t = 5, r = 0, b = 0, l = 0)),
        legend.title = element_text(size = 11), legend.text=element_text(size=11),
        legend.key.size = unit(.6, 'cm'),
        #legend.position = c(0.8, 0.2),
        legend.spacing.y = unit(.01, 'lines'),
        plot.title = element_text(size=11, hjust = 0.5, colour="black"),
        #panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        plot.margin=unit(c(0,.01,.3,0), "cm")) + 
  scale_colour_manual(values=c("#F8766D", "#00BFC4"), name="Incumbent Running?", 
                      labels=c("No","Yes"), guide = guide_legend(title.position = "top",
                                                                 reverse=TRUE)) +
  scale_fill_manual(values=c("#F8766D", "#00BFC4"), name="Incumbent Running?", 
                    labels=c("No","Yes"), guide = guide_legend(title.position = "top",
                                                               reverse=TRUE)) +
  scale_linetype_manual(values=c("longdash","solid"), name="Incumbent Running?",
                        labels=c("No","Yes"), guide = guide_legend(title.position = "top",
                                                                   reverse=TRUE))

int.s.g

dev.off()

################################################################################
## COMBINE INTERACTION GRAPHS
################################################################################

# Figure 2. Job Approval Rating × Incumbency Status Interaction (see article)

tiff("inter.tiff", units="in", width=10, height=6, res=1200)

inter <- ggarrange(int.e.g, int.s.g, ncol = 2, nrow = 1)

inter <- grid_arrange_shared_legend(int.e.g, int.s.g, ncol = 2, nrow = 1, position='bottom')

inter

dev.off()


################################################################################
## LOAD CANDIDATES' PICUTRES/IMAGES
################################################################################

dem <- readPNG("dlogo.png")
dem <- rasterGrob(dem, interpolate=TRUE)

rep <- readPNG("rlogo.png")
rep <- rasterGrob(rep, interpolate=TRUE)

cursor <- readPNG("cursor.png")
cursor <- rasterGrob(cursor, interpolate=TRUE)

biden <- readPNG("biden1.png")
biden <- rasterGrob(biden, interpolate=TRUE)

trump <- readPNG("trump1.png")
trump <- rasterGrob(trump, interpolate=TRUE)

trump2 <- readPNG("trump2.png")
trump2 <- rasterGrob(trump2, interpolate=TRUE)

gore <- readPNG("gore.png")
gore <- rasterGrob(gore, interpolate=TRUE)

bush <- readPNG("bush.png")
bush <- rasterGrob(bush, interpolate=TRUE)

kerry <- readPNG("kerry.png")
kerry <- rasterGrob(kerry, interpolate=TRUE)

mccain <- readPNG("mccain.png")
mccain <- rasterGrob(mccain, interpolate=TRUE)

obama <- readPNG("obama.png")
obama <- rasterGrob(obama, interpolate=TRUE)

romney <- readPNG("romney.png")
romney <- rasterGrob(romney, interpolate=TRUE)

clinton <- readPNG("clinton.png")
clinton <- rasterGrob(clinton, interpolate=TRUE)


################################################################################
## 2020 ELECTION FORECAST
################################################################################

# Electoral College Forecast

ev.p <- data.frame(
  gr = c("G", "G"),
  party = c("1", "2"),
  ev = c(230, 308)
)

ev.p.g <- ggplot(ev.p, aes(x = gr, y = ev)) + 
  geom_col(aes(fill = party), show.legend=FALSE, width = .1) + 
  scale_fill_manual(values=c('#E81B23', '#0044C9'))

ev.p.g <- ev.p.g + coord_flip() + theme_void() + 
  theme(plot.margin = unit(c(0,0,-1.6,0), "inches"))

ev.p.g <- ev.p.g + annotation_custom(dem, xmin=.94, xmax=1.065, ymin=95, ymax=136)
ev.p.g <- ev.p.g + annotation_custom(biden, xmin=.94, xmax=1.35, ymin=3, ymax=74)
ev.p.g <- ev.p.g + annotation_custom(rep, xmin=.94, xmax=1.065, ymin=405, ymax=444)
ev.p.g <- ev.p.g + annotation_custom(trump, xmin=.94, xmax=1.35, ymin=315, ymax=386)
ev.p.g <- ev.p.g + annotation_custom(cursor, xmin=1, xmax=1.1, ymin=263, ymax=277)
ev.p.g <- ev.p.g + geom_text(aes(x = 1, y = 175, label = "308"), color = "white", size = 5, fontface = "bold")
ev.p.g <- ev.p.g + geom_text(aes(x = 1, y = 478, label = "230"), color = "white", size = 5, fontface = "bold")
ev.p.g <- ev.p.g + geom_text(aes(x = 1.1, y = 270, label = "270"), color = "black", size = 3.5)
ev.p.g <- ev.p.g + geom_text(aes(x = 1.17, y = 145, label = "Challenger (D)"), color = "black", size = 3, fontface = "bold")
ev.p.g <- ev.p.g + geom_text(aes(x = 1.12, y = 161, label = "47th Vice President"), color = "black", size = 3)
ev.p.g <- ev.p.g + geom_text(aes(x = 1.17, y = 461, label = "Incumbent (R)"), color = "black", size = 3, fontface = "bold")
ev.p.g <- ev.p.g + geom_text(aes(x = 1.12, y = 456, label = "45th President"), color = "black", size = 3)
ev.p.g <- ev.p.g + geom_text(aes(x = 1, y = 50, label = "Biden"), color = "white", size = 4, fontface = "bold")
ev.p.g <- ev.p.g + geom_text(aes(x = 1, y = 359, label = "Trump"), color = "white", size = 4, fontface = "bold")

ev.p.g

# Electoral College Result

ev.a <- data.frame(
  gr = c("G", "G"),
  party = c("1", "2"),
  ev = c(232, 306)
)

ev.a.g <- ggplot(ev.a, aes(x = gr, y = ev)) + 
  geom_col(aes(fill = party), show.legend=FALSE, width = .1) + 
  scale_fill_manual(values=c('#E81B23', '#0044C9'))

ev.a.g <- ev.a.g + coord_flip() + theme_void() + 
  theme(plot.margin = unit(c(0,0,-1.6,0), "inches"))

ev.a.g <- ev.a.g + annotation_custom(dem, xmin=.94, xmax=1.065, ymin=95, ymax=136)
ev.a.g <- ev.a.g + annotation_custom(biden, xmin=.94, xmax=1.35, ymin=3, ymax=74)
ev.a.g <- ev.a.g + annotation_custom(rep, xmin=.94, xmax=1.065, ymin=404, ymax=443)
ev.a.g <- ev.a.g + annotation_custom(trump, xmin=.94, xmax=1.35, ymin=310, ymax=381)
ev.a.g <- ev.a.g + annotation_custom(cursor, xmin=1, xmax=1.1, ymin=263, ymax=277)
ev.a.g <- ev.a.g + geom_text(aes(x = 1, y = 175, label = "306"), color = "white", size = 5, fontface = "bold")
ev.a.g <- ev.a.g + geom_text(aes(x = 1, y = 479, label = "232"), color = "white", size = 5, fontface = "bold")
ev.a.g <- ev.a.g + geom_text(aes(x = 1.1, y = 270, label = "270"), color = "black", size = 3.5)
ev.a.g <- ev.a.g + geom_text(aes(x = 1.17, y = 145, label = "Challenger (D)"), color = "black", size = 3, fontface = "bold")
ev.a.g <- ev.a.g + geom_text(aes(x = 1.12, y = 161, label = "47th Vice President"), color = "black", size = 3)
ev.a.g <- ev.a.g + geom_text(aes(x = 1.17, y = 461, label = "Incumbent (R)"), color = "black", size = 3, fontface = "bold")
ev.a.g <- ev.a.g + geom_text(aes(x = 1.12, y = 456, label = "45th President"), color = "black", size = 3)
ev.a.g <- ev.a.g + geom_text(aes(x = 1, y = 50, label = "Biden"), color = "white", size = 4, fontface = "bold")
ev.a.g <- ev.a.g + geom_text(aes(x = 1, y = 358, label = "Trump"), color = "white", size = 4, fontface = "bold")

ev.a.g


################################################################################
## BEFORE-THE-FACT FORECASTS, 2000-2020: EXTENDED MODEL
################################################################################

# Electoral College

df.ev <- data.frame(ev = c(271, 359, 173, 384, 227, 230, 
                           266, 286, 173, 332, 227, 232),
                    forecast = c("F", "F", "F", "F", "F", "F",
                                 "A", "A", "A", "A", "A", "A"),
                    year = c("2000", "2004", "2008", "2012", "2016", "2020")
)

# Figure 4. Before-the-Fact Forecasts and Results: Electoral College, 2000–2020 (see article)

tiff("ev_g.tiff", units="in", width=7, height=8.5, res=1200)

ev.g <- df.ev %>% ggplot(aes(x=year, y=ev, fill=year)) + 
  geom_bar(stat="identity", alpha=.25, data = . %>% filter(forecast == "A"), width = 0.9, show.legend = FALSE) +
  geom_bar(stat="identity", data = . %>% filter(forecast == "F"), width = 0.6, show.legend = FALSE) +
  scale_y_continuous("Electoral Votes", limits=c(0,400), 
                     breaks=seq(0,400,100)) +
  ggtitle("   (a) Extended Model") + geom_hline(yintercept=270, linetype="dashed") +
  scale_fill_manual(values=c("#0044C9", "#E81B23","#E81B23","#0044C9","#0044C9", "#E81B23")) + theme_minimal()

ev.g <- ev.g + 
  theme(axis.title.y=element_blank(),
        axis.title.x=element_text(margin=margin(t=8, r=0, b=0, l=0), size = 14),
        axis.text.x=element_text(size = 14),
        axis.text.y=element_text(size = 14, color="black"),
        plot.title = element_text(size = 15, margin=margin(t=0, r=0, b=8, l=0))) + coord_flip()

ev.g <- ev.g + annotation_custom(gore, xmin=.6, xmax=1.4, ymin=6, ymax=60)
ev.g <- ev.g + annotation_custom(bush, xmin=1.6, xmax=2.4, ymin=6, ymax=60)
ev.g <- ev.g + annotation_custom(mccain, xmin=2.6, xmax=3.4, ymin=6, ymax=60)
ev.g <- ev.g + annotation_custom(obama, xmin=3.6, xmax=4.4, ymin=6, ymax=60)
ev.g <- ev.g + annotation_custom(clinton, xmin=4.6, xmax=5.4, ymin=6, ymax=60)
ev.g <- ev.g + annotation_custom(trump2, xmin=5.6, xmax=6.4, ymin=6, ymax=60)
ev.g <- ev.g + geom_text(aes(x = 6.45, y = 322, label = "270 to win"), color = "black", size = 4.5)
ev.g <- ev.g + geom_text(aes(x = 6.21, y = 80, label = "230"), color = "white", size = 4)
ev.g <- ev.g + geom_text(aes(x = 6.37, y = 80, label = "232"), color = "black", size = 4)
ev.g <- ev.g + geom_text(aes(x = 5.21, y = 80, label = "227"), color = "white", size = 4)
ev.g <- ev.g + geom_text(aes(x = 5.37, y = 80, label = "227"), color = "black", size = 4)
ev.g <- ev.g + geom_text(aes(x = 4.21, y = 80, label = "384"), color = "white", size = 4)
ev.g <- ev.g + geom_text(aes(x = 4.37, y = 80, label = "332"), color = "black", size = 4)
ev.g <- ev.g + geom_text(aes(x = 3.21, y = 80, label = "173"), color = "white", size = 4)
ev.g <- ev.g + geom_text(aes(x = 3.37, y = 80, label = "173"), color = "black", size = 4)
ev.g <- ev.g + geom_text(aes(x = 2.21, y = 80, label = "359"), color = "white", size = 4)
ev.g <- ev.g + geom_text(aes(x = 2.37, y = 80, label = "286"), color = "black", size = 4)
ev.g <- ev.g + geom_text(aes(x = 1.21, y = 80, label = "271"), color = "white", size = 4)
ev.g <- ev.g + geom_text(aes(x = 1.37, y = 80, label = "266"), color = "black", size = 4)

ev.g

dev.off()

################################################################################
## BEFORE-THE-FACT FORECASTS, 2000-2020: SIMPLIFIED MODEL
################################################################################

# Electoral College

df.ev2 <- data.frame(ev = c(259, 255, 111, 332, 217, 230, 
                           266, 286, 173, 332, 227, 232),
                    forecast = c("F", "F", "F", "F", "F", "F",
                                 "A", "A", "A", "A", "A", "A"),
                    year = c("2000", "2004", "2008", "2012", "2016", "2020")
)

tiff("ev_a_g.tiff", units="in", width=7, height=8.5, res=1200)

ev2.g <- df.ev2 %>% ggplot(aes(x=year, y=ev, fill=year)) + 
  geom_bar(stat="identity", alpha=.25, data = . %>% filter(forecast == "A"), width = 0.9, show.legend = FALSE) +
  geom_bar(stat="identity", data = . %>% filter(forecast == "F"), width = 0.6, show.legend = FALSE) +
  scale_y_continuous("Electoral Votes", limits=c(0,400), 
                     breaks=seq(0,400,100)) +
  ggtitle("   (b) Simplified Model") + geom_hline(yintercept=270, linetype="dashed") +
  scale_fill_manual(values=c("#0044C9", "#E81B23","#E81B23","#0044C9","#0044C9", "#E81B23")) + theme_minimal()

ev2.g <- ev2.g + 
  theme(axis.title.y=element_blank(),
        axis.title.x=element_text(margin=margin(t=8, r=0, b=0, l=0), size = 14),
        axis.text.x=element_text(size = 14),
        axis.text.y=element_text(size = 14, color="white"),
        plot.title = element_text(size = 15, margin=margin(t=0, r=0, b=8, l=0))) + coord_flip()

ev2.g <- ev2.g + annotation_custom(gore, xmin=.6, xmax=1.4, ymin=6, ymax=60)
ev2.g <- ev2.g + annotation_custom(bush, xmin=1.6, xmax=2.4, ymin=6, ymax=60)
ev2.g <- ev2.g + annotation_custom(mccain, xmin=2.6, xmax=3.4, ymin=6, ymax=60)
ev2.g <- ev2.g + annotation_custom(obama, xmin=3.6, xmax=4.4, ymin=6, ymax=60)
ev2.g <- ev2.g + annotation_custom(clinton, xmin=4.6, xmax=5.4, ymin=6, ymax=60)
ev2.g <- ev2.g + annotation_custom(trump2, xmin=5.6, xmax=6.4, ymin=6, ymax=60)
ev2.g <- ev2.g + geom_text(aes(x = 6.45, y = 322, label = "270 to win"), color = "black", size = 4.5)
ev2.g <- ev2.g + geom_text(aes(x = 6.21, y = 80, label = "230"), color = "white", size = 4)
ev2.g <- ev2.g + geom_text(aes(x = 6.37, y = 80, label = "232"), color = "black", size = 4)
ev2.g <- ev2.g + geom_text(aes(x = 5.21, y = 80, label = "217"), color = "white", size = 4)
ev2.g <- ev2.g + geom_text(aes(x = 5.37, y = 80, label = "227"), color = "black", size = 4)
ev2.g <- ev2.g + geom_text(aes(x = 4.21, y = 80, label = "332"), color = "white", size = 4)
ev2.g <- ev2.g + geom_text(aes(x = 4.37, y = 80, label = "332"), color = "black", size = 4)
ev2.g <- ev2.g + geom_text(aes(x = 3.21, y = 80, label = "111"), color = "white", size = 4)
ev2.g <- ev2.g + geom_text(aes(x = 3.37, y = 80, label = "173"), color = "black", size = 4)
ev2.g <- ev2.g + geom_text(aes(x = 2.21, y = 80, label = "255"), color = "white", size = 4)
ev2.g <- ev2.g + geom_text(aes(x = 2.37, y = 80, label = "286"), color = "black", size = 4)
ev2.g <- ev2.g + geom_text(aes(x = 1.21, y = 80, label = "259"), color = "white", size = 4)
ev2.g <- ev2.g + geom_text(aes(x = 1.37, y = 80, label = "266"), color = "black", size = 4)

ev2.g

dev.off()

# Combine graphs (Electoral College)

tiff("forecasts.tiff", units="in", width=10.5, height=7.5, res=1200)

forecasts <- ggarrange(ev.g, ev2.g, ncol = 2, nrow = 1)
forecasts

dev.off()

################################################################################
## BEFORE-THE-FACT FORECASTS, 2000-2020: EXTENDED MODEL, MULTILEVEL
################################################################################

# Electoral College

df.ev <- data.frame(ev = c(373, 345, 286, 332, 206, 356, 
                           266, 286, 173, 332, 227, 232),
                    forecast = c("F", "F", "F", "F", "F", "F",
                                 "A", "A", "A", "A", "A", "A"),
                    year = c("2000", "2004", "2008", "2012", "2016", "2020")
)

# Figure H1. Before-the-fact forecasts and results: Electoral College, 2000–2020 (see Appendix)

tiff("ev_g.tiff", units="in", width=7, height=8.5, res=1200)

ev.g <- df.ev %>% ggplot(aes(x=year, y=ev, fill=year)) + 
  geom_bar(stat="identity", alpha=.25, data = . %>% filter(forecast == "A"), width = 0.9, show.legend = FALSE) +
  geom_bar(stat="identity", data = . %>% filter(forecast == "F"), width = 0.6, show.legend = FALSE) +
  scale_y_continuous("Electoral Votes", limits=c(0,400), 
                     breaks=seq(0,400,100)) +
  ggtitle("   (a) Extended Model") + geom_hline(yintercept=270, linetype="dashed") +
  scale_fill_manual(values=c("#0044C9", "#E81B23","#E81B23","#0044C9","#0044C9", "#E81B23")) + theme_minimal()

ev.g <- ev.g + 
  theme(axis.title.y=element_blank(),
        axis.title.x=element_text(margin=margin(t=8, r=0, b=0, l=0), size = 14),
        axis.text.x=element_text(size = 14),
        axis.text.y=element_text(size = 14, color="black"),
        plot.title = element_text(size = 15, margin=margin(t=0, r=0, b=8, l=0))) + coord_flip()

ev.g <- ev.g + annotation_custom(gore, xmin=.6, xmax=1.4, ymin=6, ymax=60)
ev.g <- ev.g + annotation_custom(bush, xmin=1.6, xmax=2.4, ymin=6, ymax=60)
ev.g <- ev.g + annotation_custom(mccain, xmin=2.6, xmax=3.4, ymin=6, ymax=60)
ev.g <- ev.g + annotation_custom(obama, xmin=3.6, xmax=4.4, ymin=6, ymax=60)
ev.g <- ev.g + annotation_custom(clinton, xmin=4.6, xmax=5.4, ymin=6, ymax=60)
ev.g <- ev.g + annotation_custom(trump2, xmin=5.6, xmax=6.4, ymin=6, ymax=60)
ev.g <- ev.g + geom_text(aes(x = 6.45, y = 322, label = "270 to win"), color = "black", size = 4.5)
ev.g <- ev.g + geom_text(aes(x = 6.21, y = 80, label = "356"), color = "white", size = 4)
ev.g <- ev.g + geom_text(aes(x = 6.37, y = 80, label = "232"), color = "black", size = 4)
ev.g <- ev.g + geom_text(aes(x = 5.21, y = 80, label = "206"), color = "white", size = 4)
ev.g <- ev.g + geom_text(aes(x = 5.37, y = 80, label = "227"), color = "black", size = 4)
ev.g <- ev.g + geom_text(aes(x = 4.21, y = 80, label = "332"), color = "white", size = 4)
ev.g <- ev.g + geom_text(aes(x = 4.37, y = 80, label = "332"), color = "black", size = 4)
ev.g <- ev.g + geom_text(aes(x = 3.21, y = 80, label = "286"), color = "white", size = 4)
ev.g <- ev.g + geom_text(aes(x = 3.37, y = 80, label = "173"), color = "black", size = 4)
ev.g <- ev.g + geom_text(aes(x = 2.21, y = 80, label = "345"), color = "white", size = 4)
ev.g <- ev.g + geom_text(aes(x = 2.37, y = 80, label = "286"), color = "black", size = 4)
ev.g <- ev.g + geom_text(aes(x = 1.21, y = 80, label = "373"), color = "white", size = 4)
ev.g <- ev.g + geom_text(aes(x = 1.37, y = 80, label = "266"), color = "black", size = 4)

ev.g

dev.off()

################################################################################
## BEFORE-THE-FACT FORECASTS, 2000-2020: SIMPLIFIED MODEL, MULTILEVEL
################################################################################

# Electoral College

df.ev2 <- data.frame(ev = c(352, 251, 200, 332, 211, 212, 
                            266, 286, 173, 332, 227, 232),
                     forecast = c("F", "F", "F", "F", "F", "F",
                                  "A", "A", "A", "A", "A", "A"),
                     year = c("2000", "2004", "2008", "2012", "2016", "2020")
)

tiff("ev_a_g.tiff", units="in", width=7, height=8.5, res=1200)

ev2.g <- df.ev2 %>% ggplot(aes(x=year, y=ev, fill=year)) + 
  geom_bar(stat="identity", alpha=.25, data = . %>% filter(forecast == "A"), width = 0.9, show.legend = FALSE) +
  geom_bar(stat="identity", data = . %>% filter(forecast == "F"), width = 0.6, show.legend = FALSE) +
  scale_y_continuous("Electoral Votes", limits=c(0,400), 
                     breaks=seq(0,400,100)) +
  ggtitle("   (b) Simplified Model") + geom_hline(yintercept=270, linetype="dashed") +
  scale_fill_manual(values=c("#0044C9", "#E81B23","#E81B23","#0044C9","#0044C9", "#E81B23")) + theme_minimal()

ev2.g <- ev2.g + 
  theme(axis.title.y=element_blank(),
        axis.title.x=element_text(margin=margin(t=8, r=0, b=0, l=0), size = 14),
        axis.text.x=element_text(size = 14),
        axis.text.y=element_text(size = 14, color="white"),
        plot.title = element_text(size = 15, margin=margin(t=0, r=0, b=8, l=0))) + coord_flip()

ev2.g <- ev2.g + annotation_custom(gore, xmin=.6, xmax=1.4, ymin=6, ymax=60)
ev2.g <- ev2.g + annotation_custom(bush, xmin=1.6, xmax=2.4, ymin=6, ymax=60)
ev2.g <- ev2.g + annotation_custom(mccain, xmin=2.6, xmax=3.4, ymin=6, ymax=60)
ev2.g <- ev2.g + annotation_custom(obama, xmin=3.6, xmax=4.4, ymin=6, ymax=60)
ev2.g <- ev2.g + annotation_custom(clinton, xmin=4.6, xmax=5.4, ymin=6, ymax=60)
ev2.g <- ev2.g + annotation_custom(trump2, xmin=5.6, xmax=6.4, ymin=6, ymax=60)
ev2.g <- ev2.g + geom_text(aes(x = 6.45, y = 322, label = "270 to win"), color = "black", size = 4.5)
ev2.g <- ev2.g + geom_text(aes(x = 6.21, y = 80, label = "212"), color = "white", size = 4)
ev2.g <- ev2.g + geom_text(aes(x = 6.37, y = 80, label = "232"), color = "black", size = 4)
ev2.g <- ev2.g + geom_text(aes(x = 5.21, y = 80, label = "211"), color = "white", size = 4)
ev2.g <- ev2.g + geom_text(aes(x = 5.37, y = 80, label = "227"), color = "black", size = 4)
ev2.g <- ev2.g + geom_text(aes(x = 4.21, y = 80, label = "332"), color = "white", size = 4)
ev2.g <- ev2.g + geom_text(aes(x = 4.37, y = 80, label = "332"), color = "black", size = 4)
ev2.g <- ev2.g + geom_text(aes(x = 3.21, y = 80, label = "200"), color = "white", size = 4)
ev2.g <- ev2.g + geom_text(aes(x = 3.37, y = 80, label = "173"), color = "black", size = 4)
ev2.g <- ev2.g + geom_text(aes(x = 2.21, y = 80, label = "251"), color = "white", size = 4)
ev2.g <- ev2.g + geom_text(aes(x = 2.37, y = 80, label = "286"), color = "black", size = 4)
ev2.g <- ev2.g + geom_text(aes(x = 1.21, y = 80, label = "352"), color = "white", size = 4)
ev2.g <- ev2.g + geom_text(aes(x = 1.37, y = 80, label = "266"), color = "black", size = 4)

ev2.g

dev.off()

# Combine graphs (Electoral College)

tiff("forecasts_multilevel.tiff", units="in", width=10.5, height=7.5, res=1200)

forecasts <- ggarrange(ev.g, ev2.g, ncol = 2, nrow = 1)
forecasts

dev.off()

################################################################################
## CORRELATION BETWEEN INCUMBENT AND THIRD-PARTY VOTE
################################################################################

independents <- read_excel("independents.xlsx", sheet="Feuil1")

independents$party <- as.factor(independents$party)
independents$pinc <- as.numeric(independents$pinc)
independents$pind <- as.numeric(independents$pind)

# Figure A2. Association between incumbent and independent/third-party national vote shares, 1980–2020 (see Appendix)

tiff("scatter_ind.tiff", units="in", width=7.5, height=6, res=1200)

scatter_ind.g <- ggplot(independents, aes(x=pind, y=pinc)) + 
  geom_point(aes(color=party), size=2) +
  geom_smooth(method=lm, se=FALSE, color="grey", alpha=.7, size=.5) +
  geom_text(aes(label=year), hjust = -.3, size=3.5) +  
  scale_x_continuous("Independent and Third-Party Vote (%)", 
                     limits=c(0,20), breaks = seq(0,20,5)) +
  scale_y_continuous("Incumbent Vote (%)", 
                     limits=c(35,60), breaks = seq(35,60,5)) +
  stat_cor(method = "pearson", label.x = 16.4, label.y = 60, size = 3.5) +
  scale_color_manual("Incumbent", values=c("#0C56BC", "red3")) +
  guides(colour = guide_legend(title.position = "top")) +
  theme_bw() + theme(axis.title = element_text(size = 12.5), 
                     axis.text = element_text(size = 12),
                     axis.title.x = element_text(margin = margin(t = 5, r = 0, b = 0, l = 0)),
                     axis.title.y = element_text(margin = margin(t = 0, r = 5, b = 0, l = 0)),
                     legend.position="bottom", legend.title = element_text(hjust = 0.5, size=12),
                     legend.text=element_text(size=12))

scatter_ind.g <- scatter_ind.g + 
  #annotate("text", x = 13.65, y = 42.3, label = "(RCP polls, June)", size=3.5) +
  annotate("text", x = 9, y = 40.1, label = "Anderson: 6.6%", size=3) +
  annotate("text", x = 19.8, y = 36.5, label = "Perot: 18.9%", size=3) +
  annotate("text", x = 11, y = 48.3, label = "Perot (Reform): 8.4%", size=3)
  #annotate("text", x = 4.5, y = 49.35, label = "Nader (Green): 2.7%", size=3) +
  #annotate("text", x = 12.1, y = 41.3, label = "Kennedy: 7.6%", size=3)
  
scatter_ind.g

dev.off()


################################################################################
## LOAD FORECAST DATA
################################################################################

# Load data (extended model)

forecasts_1980_2024_e <- read_dta("forecasts_1980_2024_extended.dta")

# Load data (simplified model)

forecasts_1980_2024_s <- read_dta("forecasts_1980_2024_simplified.dta")

# Create separate dataframes per election (extended model)

forecasts_1980_e <- forecasts_1980_2024_e %>% filter(election == 1980)
forecasts_1984_e <- forecasts_1980_2024_e %>% filter(election == 1984)
forecasts_1988_e <- forecasts_1980_2024_e %>% filter(election == 1988)
forecasts_1992_e <- forecasts_1980_2024_e %>% filter(election == 1992)
forecasts_1996_e <- forecasts_1980_2024_e %>% filter(election == 1996)
forecasts_2000_e <- forecasts_1980_2024_e %>% filter(election == 2000)
forecasts_2004_e <- forecasts_1980_2024_e %>% filter(election == 2004)
forecasts_2008_e <- forecasts_1980_2024_e %>% filter(election == 2008)
forecasts_2012_e <- forecasts_1980_2024_e %>% filter(election == 2012)
forecasts_2016_e <- forecasts_1980_2024_e %>% filter(election == 2016)
forecasts_2020_e <- forecasts_1980_2024_e %>% filter(election == 2020)
forecasts_2024_e <- forecasts_1980_2024_e %>% filter(election == 2024)

# Create separate dataframes per election (simplified model)

forecasts_1980_s <- forecasts_1980_2024_s %>% filter(election == 1980)
forecasts_1984_s <- forecasts_1980_2024_s %>% filter(election == 1984)
forecasts_1988_s <- forecasts_1980_2024_s %>% filter(election == 1988)
forecasts_1992_s <- forecasts_1980_2024_s %>% filter(election == 1992)
forecasts_1996_s <- forecasts_1980_2024_s %>% filter(election == 1996)
forecasts_2000_s <- forecasts_1980_2024_s %>% filter(election == 2000)
forecasts_2004_s <- forecasts_1980_2024_s %>% filter(election == 2004)
forecasts_2008_s <- forecasts_1980_2024_s %>% filter(election == 2008)
forecasts_2012_s <- forecasts_1980_2024_s %>% filter(election == 2012)
forecasts_2016_s <- forecasts_1980_2024_s %>% filter(election == 2016)
forecasts_2020_s <- forecasts_1980_2024_s %>% filter(election == 2020)
forecasts_2024_s <- forecasts_1980_2024_s %>% filter(election == 2024)

forecasts_1980_2024_Biden <- read_dta("forecasts_1980_2024_Biden.dta")

forecasts_2024_Biden <- forecasts_1980_2024_Biden %>% filter(election == 2024)


################################################################################
## ERRORS BY STATE, 1980-2020
################################################################################

errors_e <- forecasts_1980_2024_e %>% filter(election!=2024) #remove 2024 election
errors_s <- forecasts_1980_2024_s %>% filter(election!=2024) #remove 2024 election

errors_e <- subset(errors_e, select = c(model,incumbent,election,state,abbr,err_o,mabserr_state_o,twoincv,ftwoincv_o))

View(errors_e) # data used for Figure 3 (the "mabserr_state_o" column shows the mean absolute error of out-of-sample forecasts by state which appear in Figure 3 in the article )

errors_s <- subset(errors_s, select = c(model,incumbent,election,state,abbr,err_o,mabserr_state_o,twoincv,ftwoincv_o))

errors <- rbind(errors_e,errors_s)

errors$model <- as.factor(errors$model)

# Raw error by state and election year

tiff("errors.tiff", units="in", width=8, height=8, res=1200)

errors.g <- ggplot(errors) +
  geom_point(aes(x=err_o, y=abbr, color=model, size=2)) + scale_y_discrete("", guide = guide_axis(n.dodge = 2)) +
  geom_vline(xintercept = 0) + facet_wrap(~ election, nrow=2) + 
  scale_x_continuous("Error", limits=c(-25,25), breaks = seq(-25,25,5)) + 
  theme_bw()

errors.g

dev.off()

# Mean absolute error by state (see Figure 3 in the article)

tiff("mae_state.tiff", units="in", width=5, height=7, res=1200)

mae.state.g <- errors_e %>%
  arrange(mabserr_state_o) %>%
  mutate(abbr = forcats::fct_inorder(abbr)) %>% ggplot() +
  geom_point(aes(x=mabserr_state_o, y=abbr), size=1) + 
  scale_y_discrete("", guide = guide_axis(n.dodge = 2)) +
  geom_vline(xintercept=3.58, linetype="dashed", color="black", size=.3) +
  #geom_vline(xintercept=3.47, linetype="dashed", color="#56B4E9", size=.3) +
  #scale_color_manual(values=c("#E69F00", "#56B4E9")) +
  scale_x_continuous("Mean Absolute Error", limits=c(0,8), breaks = seq(0,8,2)) +
  theme_bw() + theme(axis.title.x = element_text(size=9, 
               margin = margin(t = 3, r = 0, b = 0, l = 0)),
               axis.text.x = element_text(size=7),
               axis.text.y = element_text(size=7))

mae.state.g <- mae.state.g + geom_text(aes(x = 2.4, y = 50.3, label = "Overall MAE"), color = "black", size = 2.5)

mae.state.g

dev.off()

# Mean absolute error by state, map

tiff("mae_map_e.tiff", units="in", width=5, height=5, res=1200)

mae.map.e <- plot_usmap(regions = "state", data = errors_e, 
                      values = "mabserr_state_o", labels = FALSE, color="#d24e01") + 
                      scale_fill_gradientn("Mean Absolute Error", colours = RColorBrewer::brewer.pal(9, "Oranges"),
                      na.value = "grey90",
                      guide = guide_colourbar(barwidth = 15, 
                                               barheight = 0.4,
                                               #put legend title on top of legend 
                                               title.position = "top")) +
                      theme(plot.title = element_blank()) + 
                      theme(legend.position = "bottom", legend.text = element_text(size = 8, 
                      margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
                      "lines")) + theme(legend.key.size = unit(0.4, "cm")) + 
                      theme(legend.position = "bottom"); mae.map.e

dev.off()

tiff("mae_map_s.tiff", units="in", width=5, height=5, res=1200)

mae.map.s <- plot_usmap(regions = "state", data = errors_s, 
                        values = "mabserr_state_o", labels = FALSE, color="#d24e01") + 
  scale_fill_gradientn("Mean Absolute Error", colours = RColorBrewer::brewer.pal(9, "Oranges"),
                       na.value = "grey90",
                       guide = guide_colourbar(barwidth = 15, 
                                               barheight = 0.4,
                                               #put legend title on top of legend 
                                               title.position = "top")) +
  theme(plot.title = element_blank()) + 
  theme(legend.position = "bottom", legend.text = element_text(size = 8, 
  margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
                                                                                                                                 "lines")) + theme(legend.key.size = unit(0.4, "cm")) + 
  theme(legend.position = "bottom"); mae.map.s

dev.off()

# Combine graphs (MAE by state and MAE map)

# Figure 3. Mean Absolute Error by State From Out-of-Sample Forecasts (Extended Model), 1980–2020 (see article)

tiff("mae_1980_2020.tiff", units="in", width=7, height=4, res=1200)

mae.1980.2020 <- ggarrange(mae.state.g, mae.map.e, ncol = 2); mae.1980.2020

dev.off()

# Mean absolute error by state and election year

errors.us.e <- forecasts_1980_2024_e %>% filter(election!=2024) #remove 2024 election

errors.us.e <- subset(errors.us.e, select = c(election, state, abbr, abserr_o))

errors.us.e <- errors.us.e %>% 
  sjlabelled::remove_all_labels() %>% 
  tibble()

errors.us.e <- errors.us.e %>% 
  rename(year = election,
         state.name = state,
         state = abbr,
         n = abserr_o)

errors.us.e$state.name[errors_e$state == "District of Columbia"] <- "Washington DC"

bl = sapply(1:37, function(n) paste(rep(" ",n),collapse="")) #create unique blank strip labels for empty facets

errors.us.e$state.reorder <- factor(errors.us.e$state.name,
                                  levels = c(bl[1:10], "Maine",
                                             bl[11:19], "Vermont", "New Hampshire",
                                             "Washington", "Idaho", "Montana", "North Dakota", "Minnesota", "Illinois", "Wisconsin", "Michigan", "New York", "Massachusetts", "Rhode Island",
                                             "Oregon", "Nevada", "Wyoming", "South Dakota", "Iowa", "Indiana", "Ohio", "Pennsylvania", "New Jersey", "Connecticut", bl[20],
                                             "California", "Utah", "Colorado", "Nebraska", "Missouri", "Kentucky", "West Virginia", "Virginia", "Maryland", "Washington DC", bl[21],
                                             bl[22], "Arizona", "New Mexico", "Kansas", "Arkansas", "Tennessee", "North Carolina", "South Carolina", "Delaware", bl[23:24],
                                             bl[25:27], "Oklahoma", "Louisiana", "Mississippi", "Alabama", "Georgia", bl[28:29],
                                             bl[30], "Alaska", "Hawaii",
                                             bl[31], "Texas", bl[32:35], "Florida"))

# Figure E1. Out-of-sample mean absolute errors by election year and state, 1980–2020 – Extended model (see Appendix)

tiff("errors_us_e.tiff", units="in", width=12, height=9, res=1200)

errors.us.e.g <- ggplot(subset(errors.us.e, is.na(state.name)==F), aes(y = n, x = year)) +
  geom_line() + facet_wrap( ~ state.reorder, ncol = 11, drop = F, strip.position="bottom") +
  scale_y_continuous("Mean Absolute Error", 
                     limits=c(0,30), breaks = seq(0,30,10)) +
  theme_classic() +
  theme(axis.text.y = element_text(size=12),
        axis.text.x = element_blank(),
        axis.title.x=element_text(margin=margin(t=8, r=0, b=0, l=0), size = 14),
        axis.title.y=element_text(margin=margin(t=0, r=8, b=0, l=0), size = 14),
        strip.background=element_blank(),
        axis.line=element_blank(),
        axis.ticks=element_blank()) +
  ylab("Mean Absolute Error") +
  xlab("Election (1980–2020)")

errors.us.e.g

dev.off()

errors.us.s <- forecasts_1980_2024_s %>% filter(election!=2024) #remove 2024 election

errors.us.s <- subset(errors.us.s, select = c(election, state, abbr, abserr_o))

errors.us.s <- errors.us.s %>% 
  sjlabelled::remove_all_labels() %>% 
  tibble()

errors.us.s <- errors.us.s %>% 
  rename(year = election,
         state.name = state,
         state = abbr,
         n = abserr_o)

errors.us.s$state.name[errors_s$state == "District of Columbia"] <- "Washington DC"

bl = sapply(1:37, function(n) paste(rep(" ",n),collapse="")) #create unique blank strip labels for empty facets

errors.us.s$state.reorder <- factor(errors.us.s$state.name,
                                    levels = c(bl[1:10], "Maine",
                                               bl[11:19], "Vermont", "New Hampshire",
                                               "Washington", "Idaho", "Montana", "North Dakota", "Minnesota", "Illinois", "Wisconsin", "Michigan", "New York", "Massachusetts", "Rhode Island",
                                               "Oregon", "Nevada", "Wyoming", "South Dakota", "Iowa", "Indiana", "Ohio", "Pennsylvania", "New Jersey", "Connecticut", bl[20],
                                               "California", "Utah", "Colorado", "Nebraska", "Missouri", "Kentucky", "West Virginia", "Virginia", "Maryland", "Washington DC", bl[21],
                                               bl[22], "Arizona", "New Mexico", "Kansas", "Arkansas", "Tennessee", "North Carolina", "South Carolina", "Delaware", bl[23:24],
                                               bl[25:27], "Oklahoma", "Louisiana", "Mississippi", "Alabama", "Georgia", bl[28:29],
                                               bl[30], "Alaska", "Hawaii",
                                               bl[31], "Texas", bl[32:35], "Florida"))

# Figure E2. Out-of-sample mean absolute errors by election year and state, 1980–2020 – Simplified model (see Appendix)

tiff("errors_us_s.tiff", units="in", width=12, height=9, res=1200)

errors.us.s.g <- ggplot(subset(errors.us.s, is.na(state.name)==F), aes(y = n, x = year)) +
  geom_line() + facet_wrap( ~ state.reorder, ncol = 11, drop = F, strip.position="bottom") +
  scale_y_continuous("Mean Absolute Error", 
                     limits=c(0,30), breaks = seq(0,30,10)) +
  theme_classic() +
  theme(axis.text.y = element_text(size=12),
        axis.text.x = element_blank(),
        axis.title.x=element_text(margin=margin(t=8, r=0, b=0, l=0), size = 14),
        axis.title.y=element_text(margin=margin(t=0, r=8, b=0, l=0), size = 14),
        strip.background=element_blank(),
        axis.line=element_blank(),
        axis.ticks=element_blank()) +
  ylab("Mean Absolute Error") +
  xlab("Election (1980–2020)")

errors.us.s.g

dev.off()

# Scatter plot: predicted vs observed values

tiff("scatter.tiff", units="in", width=8, height=7, res=1200)

scatter.g <- ggplot(errors_e, aes(x=twoincv, y=ftwoincv_o)) +
  geom_point(size=2, alpha = 0.5) + facet_wrap(~ election, nrow=3) +
  geom_smooth(method=lm, se=FALSE, color="darkred") + 
  scale_x_continuous("Observed", limits=c(0,100), breaks = seq(0,100,20)) + 
  scale_y_continuous("Predicted", limits=c(0,100), breaks = seq(0,100,20)) +
  theme_bw()

scatter.g

dev.off()

# Error distribution by party

tiff("distribution_ridges_errors.tiff", units="in", width=5, height=5, res=1200)

distribution.ridges.errors.g <- ggplot(errors_e, aes(x = err_o, y = incumbent, fill = incumbent)) +
  geom_density_ridges(alpha = 0.7) +
  theme_ridges() +
  geom_vline(xintercept=0, linetype="dashed", color="black", size=.3) +
  coord_cartesian(clip = "off") +
  theme_bw() +
  theme(legend.position = "none",
        axis.title = element_text(size = 12), 
        axis.text = element_text(size = 12),
        axis.title.x = element_text(margin = margin(t = 5, r = 0, b = 0, l = 0))) + 
  scale_x_continuous("Error") +
  scale_y_discrete("", expand = c(-0.02, 0.3),
                   labels=c("DEM"="Democrat","REP"="Republican")) + 
  scale_fill_manual(values = c("#0C56BC", "red3"))

distribution.ridges.errors.g

dev.off()

################################################################################
## 2000 ELECTION: EXTENDED MODEL
################################################################################

# Figure F1. Predicted and actual outcomes, 2000 presidential election (see Appendix)

forecast.map.e.2000 <- plot_usmap(regions = "state", data = forecasts_2000_e, 
values = "fstatewinner_b", labels = FALSE) + scale_fill_manual(name = "
                   ", values = c("#0C56BC", "red3")) + 
theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
theme(legend.position = "bottom", legend.text = element_text(size = 8, 
margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
"lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(a) Predicted winner in each state"); 
forecast.map.e.2000

outcome.map.e.2000 <- plot_usmap(regions = "state", data = forecasts_2000_e, 
values = "statewinner", labels = FALSE) + scale_fill_manual(name = "
                    ", values = c("#0C56BC", "red3")) + 
theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
theme(legend.position = "bottom", legend.text = element_text(size = 8, 
margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
"lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(b) Actual winner in each state"); 
outcome.map.e.2000

correct.map.e.2000 <- plot_usmap(regions = "state", data = forecasts_2000_e, 
values = "statecorrect_b", labels = FALSE) + scale_fill_manual(name = "
                             ", values = c("#00b300", "#e5e5e5"), limits = c("Correct", "Incorrect")) +  
theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
theme(legend.position = "bottom", legend.text = element_text(size = 8, 
margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
"lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(c) Correct and incorrect forecasts"); 
correct.map.e.2000

layout_matrix <- matrix(c(1, 1, 2, 2, 4, 3, 3, 4), nrow = 2, byrow = TRUE)
map.e.2000 <- grid.arrange(forecast.map.e.2000, outcome.map.e.2000, correct.map.e.2000, 
layout_matrix = layout_matrix); map.e.2000

dev.copy(pdf, "map_e_2000.pdf")
dev.off()


################################################################################
## 2004 ELECTION: EXTENDED MODEL
################################################################################

# Figure F2. Predicted and actual outcomes, 2004 presidential election (see Appendix)

forecast.map.e.2004 <- plot_usmap(regions = "state", data = forecasts_2004_e, 
values = "fstatewinner_b", labels = FALSE) + scale_fill_manual(name = "
                  ", values = c("red3", "#0C56BC")) + 
theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
theme(legend.position = "bottom", legend.text = element_text(size = 8, 
margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
"lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(a) Predicted winner in each state"); 
forecast.map.e.2004

outcome.map.e.2004 <- plot_usmap(regions = "state", data = forecasts_2004_e, 
values = "statewinner", labels = FALSE) + scale_fill_manual(name = "
                   ", values = c("red3", "#0C56BC")) + 
theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
theme(legend.position = "bottom", legend.text = element_text(size = 8, 
margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
"lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(b) Actual winner in each state"); 
outcome.map.e.2004

correct.map.e.2004 <- plot_usmap(regions = "state", data = forecasts_2004_e, 
values = "statecorrect_b", labels = FALSE) + scale_fill_manual(name = "
                             ", values = c("#00b300", "#e5e5e5"), limits = c("Correct", "Incorrect")) +  
theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
theme(legend.position = "bottom", legend.text = element_text(size = 8, 
margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
"lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(c) Correct and incorrect forecasts"); 
correct.map.e.2004

layout_matrix <- matrix(c(1, 1, 2, 2, 4, 3, 3, 4), nrow = 2, byrow = TRUE)
map.e.2004 <- grid.arrange(forecast.map.e.2004, outcome.map.e.2004, correct.map.e.2004, 
layout_matrix = layout_matrix); map.e.2004

dev.copy(pdf, "map_e_2004.pdf")
dev.off()


################################################################################
## 2008 ELECTION: EXTENDED MODEL
################################################################################

# Figure F3. Predicted and actual outcomes, 2008 presidential election (see Appendix)

forecast.map.e.2008 <- plot_usmap(regions = "state", data = forecasts_2008_e, 
values = "fstatewinner_b", labels = FALSE) + scale_fill_manual(name = "
                  ", values = c("#0C56BC", "red3")) + 
theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
theme(legend.position = "bottom", legend.text = element_text(size = 8, 
margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
"lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(a) Predicted winner in each state"); 
forecast.map.e.2008

outcome.map.e.2008 <- plot_usmap(regions = "state", data = forecasts_2008_e, 
values = "statewinner", labels = FALSE) + scale_fill_manual(name = "
                   ", values = c("#0C56BC", "red3")) + 
theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
theme(legend.position = "bottom", legend.text = element_text(size = 8, 
margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
"lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(b) Actual winner in each state"); 
outcome.map.e.2008

correct.map.e.2008 <- plot_usmap(regions = "state", data = forecasts_2008_e, 
values = "statecorrect_b", labels = FALSE) + scale_fill_manual(name = "
                             ", values = c("#00b300", "#e5e5e5"), limits = c("Correct", "Incorrect")) +  
theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
theme(legend.position = "bottom", legend.text = element_text(size = 8, 
margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
"lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(c) Correct and incorrect forecasts"); 
correct.map.e.2008

layout_matrix <- matrix(c(1, 1, 2, 2, 4, 3, 3, 4), nrow = 2, byrow = TRUE)
map.e.2008 <- grid.arrange(forecast.map.e.2008, outcome.map.e.2008, correct.map.e.2008, 
layout_matrix = layout_matrix); map.e.2008

dev.copy(pdf, "map_e_2008.pdf")
dev.off()


################################################################################
## 2012 ELECTION: EXTENDED MODEL
################################################################################

# Figure F4. Predicted and actual outcomes, 2012 presidential election (see Appendix)

forecast.map.e.2012 <- plot_usmap(regions = "state", data = forecasts_2012_e, 
values = "fstatewinner_b", labels = FALSE) + scale_fill_manual(name = "
                 ", values = c("#0C56BC", "red3")) + 
theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
theme(legend.position = "bottom", legend.text = element_text(size = 8, 
margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
"lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(a) Predicted winner in each state"); 
forecast.map.e.2012

outcome.map.e.2012 <- plot_usmap(regions = "state", data = forecasts_2012_e, 
values = "statewinner", labels = FALSE) + scale_fill_manual(name = "
                   ", values = c("#0C56BC", "red3")) + 
theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
theme(legend.position = "bottom", legend.text = element_text(size = 8, 
margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
"lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(b) Actual winner in each state"); 
outcome.map.e.2012

correct.map.e.2012 <- plot_usmap(regions = "state", data = forecasts_2012_e, 
values = "statecorrect_b", labels = FALSE) + scale_fill_manual(name = "
                             ", values = c("#00b300", "#e5e5e5"), limits = c("Correct", "Incorrect")) +  
theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
theme(legend.position = "bottom", legend.text = element_text(size = 8, 
margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
"lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(c) Correct and incorrect forecasts"); 
correct.map.e.2012

layout_matrix <- matrix(c(1, 1, 2, 2, 4, 3, 3, 4), nrow = 2, byrow = TRUE)
map.e.2012 <- grid.arrange(forecast.map.e.2012, outcome.map.e.2012, correct.map.e.2012, 
layout_matrix = layout_matrix); map.e.2012

dev.copy(pdf, "map_e_2012.pdf")
dev.off()


################################################################################
## 2016 ELECTION: EXTENDED MODEL
################################################################################

# Figure F5. Predicted and actual outcomes, 2016 presidential election (see Appendix)

forecast.map.e.2016 <- plot_usmap(regions = "state", data = forecasts_2016_e, 
values = "fstatewinner_b", labels = FALSE) + scale_fill_manual(name = "
                  ", values = c("red3", "#0C56BC")) + 
theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
theme(legend.position = "bottom", legend.text = element_text(size = 8, 
margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
"lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(a) Predicted winner in each state"); 
forecast.map.e.2016

outcome.map.e.2016 <- plot_usmap(regions = "state", data = forecasts_2016_e, 
values = "statewinner", labels = FALSE) + scale_fill_manual(name = "
                   ", values = c("red3", "#0C56BC")) + 
theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
theme(legend.position = "bottom", legend.text = element_text(size = 8, 
margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
"lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(b) Actual winner in each state"); 
outcome.map.e.2016

correct.map.e.2016 <- plot_usmap(regions = "state", data = forecasts_2016_e, 
values = "statecorrect_b", labels = FALSE) + scale_fill_manual(name = "
                             ", values = c("#00b300", "#e5e5e5"), limits = c("Correct", "Incorrect")) +  
theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
theme(legend.position = "bottom", legend.text = element_text(size = 8, 
margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
"lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(c) Correct and incorrect forecasts"); 
correct.map.e.2016

layout_matrix <- matrix(c(1, 1, 2, 2, 4, 3, 3, 4), nrow = 2, byrow = TRUE)
map.e.2016 <- grid.arrange(forecast.map.e.2016, outcome.map.e.2016, correct.map.e.2016, 
layout_matrix = layout_matrix); map.e.2016

dev.copy(pdf, "map_e_2016.pdf")
dev.off()


################################################################################
## 2020 ELECTION #1: EXTENDED MODEL
################################################################################

# Figure F6. Predicted and actual outcomes, 2020 presidential election (see Appendix)

forecast.map.e.2020 <- plot_usmap(regions = "state", data = forecasts_2020_e, 
                                values = "fstatewinner_b", labels = FALSE) + scale_fill_manual(name = "
                  ", values = c("red3", "#0C56BC")) + 
theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
theme(legend.position = "bottom", legend.text = element_text(size = 8, 
margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
"lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(a) Predicted winner in each state"); 
forecast.map.e.2020

outcome.map.e.2020 <- plot_usmap(regions = "state", data = forecasts_2020_e, 
                               values = "statewinner", labels = FALSE) + scale_fill_manual(name = "
                   ", values = c("red3", "#0C56BC")) + 
theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
theme(legend.position = "bottom", legend.text = element_text(size = 8, 
margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
"lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(b) Actual winner in each state"); 
outcome.map.e.2020

correct.map.e.2020 <- plot_usmap(regions = "state", data = forecasts_2020_e, 
                               values = "statecorrect_b", labels = FALSE) + scale_fill_manual(name = "
                             ", values = c("#00b300", "#e5e5e5"), limits = c("Correct", "Incorrect")) +  
theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
theme(legend.position = "bottom", legend.text = element_text(size = 8, 
margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
"lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(c) Correct and incorrect forecasts"); 
correct.map.e.2020

layout_matrix <- matrix(c(1, 1, 2, 2, 4, 3, 3, 4), nrow = 2, byrow = TRUE)
map.e.2020 <- grid.arrange(forecast.map.e.2020, outcome.map.e.2020, correct.map.e.2020, 
                         layout_matrix = layout_matrix); map.e.2020

dev.copy(pdf, "map_e_2020.pdf")
dev.off()

################################################################################
## 2024 ELECTION: EXTENDED MODEL
################################################################################

forecast.map.e.2024 <- plot_usmap(regions = "state", data = forecasts_2024_e, 
                                  values = "fstatewinner_b", labels = FALSE) + scale_fill_manual(name = "
                  ", values = c("red3", "#0C56BC")) + 
  theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
  theme(legend.position = "bottom", legend.text = element_text(size = 8, 
                                                               margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
                                                                                                                                 "lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(a) Predicted winner in each state"); 
forecast.map.e.2024

outcome.map.e.2024 <- plot_usmap(regions = "state", data = forecasts_2024_e, 
                                 values = "statewinner", labels = FALSE) + scale_fill_manual(name = "
                   ", values = c("red3", "#0C56BC")) + 
  theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
  theme(legend.position = "bottom", legend.text = element_text(size = 8, 
                                                               margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
                                                                                                                                 "lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(b) Actual winner in each state"); 
outcome.map.e.2024

correct.map.e.2024 <- plot_usmap(regions = "state", data = forecasts_2024_e, 
                                 values = "statecorrect_b", labels = FALSE) + scale_fill_manual(name = "
                             ", values = c("#00b300", "#e5e5e5"), limits = c("Correct", "Incorrect")) +  
  theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
  theme(legend.position = "bottom", legend.text = element_text(size = 8, 
                                                               margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
                                                                                                                                 "lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(c) Correct and incorrect forecasts"); 
correct.map.e.2024

layout_matrix <- matrix(c(1, 1, 2, 2, 4, 3, 3, 4), nrow = 2, byrow = TRUE)
map.e.2024 <- grid.arrange(forecast.map.e.2024, outcome.map.e.2024, correct.map.e.2024, 
                           layout_matrix = layout_matrix); map.e.2024

dev.copy(pdf, "map_e_2024.pdf")
dev.off()

################################################################################
## 2000 ELECTION: SIMPLIFIED MODEL
################################################################################

# Figure G1. Predicted and actual outcomes, 2000 presidential election (see Appendix)

forecast.map.s.2000 <- plot_usmap(regions = "state", data = forecasts_2000_s, 
                                  values = "fstatewinner_b", labels = FALSE) + scale_fill_manual(name = "
                   ", values = c("#0C56BC", "red3")) + 
  theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
  theme(legend.position = "bottom", legend.text = element_text(size = 8, 
                                                               margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
                                                                                                                                 "lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(a) Predicted winner in each state"); 
forecast.map.s.2000

outcome.map.s.2000 <- plot_usmap(regions = "state", data = forecasts_2000_s, 
                                 values = "statewinner", labels = FALSE) + scale_fill_manual(name = "
                    ", values = c("#0C56BC", "red3")) + 
  theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
  theme(legend.position = "bottom", legend.text = element_text(size = 8, 
                                                               margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
                                                                                                                                 "lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(b) Actual winner in each state"); 
outcome.map.s.2000

correct.map.s.2000 <- plot_usmap(regions = "state", data = forecasts_2000_s, 
                                 values = "statecorrect_b", labels = FALSE) + scale_fill_manual(name = "
                             ", values = c("#00b300", "#e5e5e5"), limits = c("Correct", "Incorrect")) +  
  theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
  theme(legend.position = "bottom", legend.text = element_text(size = 8, 
                                                               margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
                                                                                                                                 "lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(c) Correct and incorrect forecasts"); 
correct.map.s.2000

layout_matrix <- matrix(c(1, 1, 2, 2, 4, 3, 3, 4), nrow = 2, byrow = TRUE)
map.s.2000 <- grid.arrange(forecast.map.s.2000, outcome.map.s.2000, correct.map.s.2000, 
                           layout_matrix = layout_matrix); map.s.2000

dev.copy(pdf, "map_s_2000.pdf")
dev.off()


################################################################################
## 2004 ELECTION: SIMPLIFIED MODEL
################################################################################

# Figure G2. Predicted and actual outcomes, 2004 presidential election (see Appendix)

forecast.map.s.2004 <- plot_usmap(regions = "state", data = forecasts_2004_s, 
                                  values = "fstatewinner_b", labels = FALSE) + scale_fill_manual(name = "
                  ", values = c("red3", "#0C56BC")) + 
  theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
  theme(legend.position = "bottom", legend.text = element_text(size = 8, 
                                                               margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
                                                                                                                                 "lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(a) Predicted winner in each state"); 
forecast.map.s.2004

outcome.map.s.2004 <- plot_usmap(regions = "state", data = forecasts_2004_s, 
                                 values = "statewinner", labels = FALSE) + scale_fill_manual(name = "
                   ", values = c("red3", "#0C56BC")) + 
  theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
  theme(legend.position = "bottom", legend.text = element_text(size = 8, 
                                                               margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
                                                                                                                                 "lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(b) Actual winner in each state"); 
outcome.map.s.2004

correct.map.s.2004 <- plot_usmap(regions = "state", data = forecasts_2004_s, 
                                 values = "statecorrect_b", labels = FALSE) + scale_fill_manual(name = "
                             ", values = c("#00b300", "#e5e5e5"), limits = c("Correct", "Incorrect")) +  
  theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
  theme(legend.position = "bottom", legend.text = element_text(size = 8, 
                                                               margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
                                                                                                                                 "lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(c) Correct and incorrect forecasts"); 
correct.map.s.2004

layout_matrix <- matrix(c(1, 1, 2, 2, 4, 3, 3, 4), nrow = 2, byrow = TRUE)
map.s.2004 <- grid.arrange(forecast.map.s.2004, outcome.map.s.2004, correct.map.s.2004, 
                           layout_matrix = layout_matrix); map.s.2004

dev.copy(pdf, "map_s_2004.pdf")
dev.off()


################################################################################
## 2008 ELECTION: SIMPLIFIED MODEL
################################################################################

# Figure G3. Predicted and actual outcomes, 2008 presidential election (see Appendix)

forecast.map.s.2008 <- plot_usmap(regions = "state", data = forecasts_2008_s, 
                                  values = "fstatewinner_b", labels = FALSE) + scale_fill_manual(name = "
                  ", values = c("#0C56BC", "red3")) + 
  theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
  theme(legend.position = "bottom", legend.text = element_text(size = 8, 
                                                               margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
                                                                                                                                 "lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(a) Predicted winner in each state"); 
forecast.map.s.2008

outcome.map.s.2008 <- plot_usmap(regions = "state", data = forecasts_2008_s, 
                                 values = "statewinner", labels = FALSE) + scale_fill_manual(name = "
                   ", values = c("#0C56BC", "red3")) + 
  theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
  theme(legend.position = "bottom", legend.text = element_text(size = 8, 
                                                               margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
                                                                                                                                 "lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(b) Actual winner in each state"); 
outcome.map.s.2008

correct.map.s.2008 <- plot_usmap(regions = "state", data = forecasts_2008_s, 
                                 values = "statecorrect_b", labels = FALSE) + scale_fill_manual(name = "
                             ", values = c("#00b300", "#e5e5e5"), limits = c("Correct", "Incorrect")) +  
  theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
  theme(legend.position = "bottom", legend.text = element_text(size = 8, 
                                                               margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
                                                                                                                                 "lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(c) Correct and incorrect forecasts"); 
correct.map.s.2008

layout_matrix <- matrix(c(1, 1, 2, 2, 4, 3, 3, 4), nrow = 2, byrow = TRUE)
map.s.2008 <- grid.arrange(forecast.map.s.2008, outcome.map.s.2008, correct.map.s.2008, 
                           layout_matrix = layout_matrix); map.s.2008

dev.copy(pdf, "map_s_2008.pdf")
dev.off()


################################################################################
## 2012 ELECTION: SIMPLIFIED MODEL
################################################################################

# Figure G4. Predicted and actual outcomes, 2012 presidential election (see Appendix)

forecast.map.s.2012 <- plot_usmap(regions = "state", data = forecasts_2012_s, 
                                  values = "fstatewinner_b", labels = FALSE) + scale_fill_manual(name = "
                 ", values = c("#0C56BC", "red3")) + 
  theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
  theme(legend.position = "bottom", legend.text = element_text(size = 8, 
                                                               margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
                                                                                                                                 "lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(a) Predicted winner in each state"); 
forecast.map.s.2012

outcome.map.s.2012 <- plot_usmap(regions = "state", data = forecasts_2012_s, 
                                 values = "statewinner", labels = FALSE) + scale_fill_manual(name = "
                   ", values = c("#0C56BC", "red3")) + 
  theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
  theme(legend.position = "bottom", legend.text = element_text(size = 8, 
                                                               margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
                                                                                                                                 "lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(b) Actual winner in each state"); 
outcome.map.s.2012

correct.map.s.2012 <- plot_usmap(regions = "state", data = forecasts_2012_s, 
                                 values = "statecorrect_b", labels = FALSE) + scale_fill_manual(name = "
                             ", values = c("#00b300", "#e5e5e5"), limits = c("Correct", "Incorrect")) +  
  theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
  theme(legend.position = "bottom", legend.text = element_text(size = 8, 
                                                               margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
                                                                                                                                 "lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(c) Correct and incorrect forecasts"); 
correct.map.s.2012

layout_matrix <- matrix(c(1, 1, 2, 2, 4, 3, 3, 4), nrow = 2, byrow = TRUE)
map.s.2012 <- grid.arrange(forecast.map.s.2012, outcome.map.s.2012, correct.map.s.2012, 
                           layout_matrix = layout_matrix); map.s.2012

dev.copy(pdf, "map_s_2012.pdf")
dev.off()


################################################################################
## 2016 ELECTION: SIMPLIFIED MODEL
################################################################################

# Figure G5. Predicted and actual outcomes, 2016 presidential election (see Appendix)

forecast.map.s.2016 <- plot_usmap(regions = "state", data = forecasts_2016_s, 
                                  values = "fstatewinner_b", labels = FALSE) + scale_fill_manual(name = "
                  ", values = c("red3", "#0C56BC")) + 
  theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
  theme(legend.position = "bottom", legend.text = element_text(size = 8, 
                                                               margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
                                                                                                                                 "lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(a) Predicted winner in each state"); 
forecast.map.s.2016

outcome.map.s.2016 <- plot_usmap(regions = "state", data = forecasts_2016_s, 
                                 values = "statewinner", labels = FALSE) + scale_fill_manual(name = "
                   ", values = c("red3", "#0C56BC")) + 
  theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
  theme(legend.position = "bottom", legend.text = element_text(size = 8, 
                                                               margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
                                                                                                                                 "lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(b) Actual winner in each state"); 
outcome.map.s.2016

correct.map.s.2016 <- plot_usmap(regions = "state", data = forecasts_2016_s, 
                                 values = "statecorrect_b", labels = FALSE) + scale_fill_manual(name = "
                             ", values = c("#00b300", "#e5e5e5"), limits = c("Correct", "Incorrect")) +  
  theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
  theme(legend.position = "bottom", legend.text = element_text(size = 8, 
                                                               margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
                                                                                                                                 "lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(c) Correct and incorrect forecasts"); 
correct.map.s.2016

layout_matrix <- matrix(c(1, 1, 2, 2, 4, 3, 3, 4), nrow = 2, byrow = TRUE)
map.s.2016 <- grid.arrange(forecast.map.s.2016, outcome.map.s.2016, correct.map.s.2016, 
                           layout_matrix = layout_matrix); map.s.2016

dev.copy(pdf, "map_s_2016.pdf")
dev.off()


################################################################################
## 2020 ELECTION: SIMPLIFIED MODEL
################################################################################

# Figure G6. Predicted and actual outcomes, 2020 presidential election (see Appendix)

forecast.map.s.2020 <- plot_usmap(regions = "state", data = forecasts_2020_s, 
                                  values = "fstatewinner_b", labels = FALSE) + scale_fill_manual(name = "
                  ", values = c("red3", "#0C56BC")) + 
  theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
  theme(legend.position = "bottom", legend.text = element_text(size = 8, 
                                                               margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
                                                                                                                                 "lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(a) Predicted winner in each state"); 
forecast.map.s.2020

outcome.map.s.2020 <- plot_usmap(regions = "state", data = forecasts_2020_s, 
                                 values = "statewinner", labels = FALSE) + scale_fill_manual(name = "
                   ", values = c("red3", "#0C56BC")) + 
  theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
  theme(legend.position = "bottom", legend.text = element_text(size = 8, 
                                                               margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
                                                                                                                                 "lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(b) Actual winner in each state"); 
outcome.map.s.2020

correct.map.s.2020 <- plot_usmap(regions = "state", data = forecasts_2020_s, 
                                 values = "statecorrect_b", labels = FALSE) + scale_fill_manual(name = "
                             ", values = c("#00b300", "#e5e5e5"), limits = c("Correct", "Incorrect")) +  
  theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
  theme(legend.position = "bottom", legend.text = element_text(size = 8, 
                                                               margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
                                                                                                                                 "lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(c) Correct and incorrect forecasts"); 
correct.map.s.2020

layout_matrix <- matrix(c(1, 1, 2, 2, 4, 3, 3, 4), nrow = 2, byrow = TRUE)
map.s.2020 <- grid.arrange(forecast.map.s.2020, outcome.map.s.2020, correct.map.s.2020, 
                           layout_matrix = layout_matrix); map.s.2020

dev.copy(pdf, "map_s_2020.pdf")
dev.off()


################################################################################
## 2020 ELECTION #2: EXTENDED MODEL
################################################################################

# Figure 1. U.S. Presidential Election Forecasts and Results, 2020 (see article)

tiff("2020_map.tiff", units="in", width=8, height=8, res=1200)

forecast.map.e.2020 <- plot_usmap(regions = "state", data = forecasts_2020_e, 
values = "fstatewinner_b", labels = FALSE) + scale_fill_manual(name = "
                  ", values = c("red3", "#0C56BC")) + 
theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
theme(legend.position = "bottom", legend.text = element_text(size = 8, 
margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
                                                                                                                                 "lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(a) Predicted winner in each state"); 
forecast.map.e.2020

outcome.map.e.2020 <- plot_usmap(regions = "state", data = forecasts_2020_e, 
values = "statewinner", labels = FALSE) + scale_fill_manual(name = "
                   ", values = c("red3", "#0C56BC")) + 
theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
theme(legend.position = "bottom", legend.text = element_text(size = 8, 
margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
                                                                                                                                 "lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(b) Actual winner in each state"); 
outcome.map.e.2020

correct.map.e.2020 <- plot_usmap(regions = "state", data = forecasts_2020_e, 
values = "statecorrect_b", labels = FALSE) + scale_fill_manual(name = "
                             ", values = c("#00b300", "#e5e5e5"), limits = c("Correct", "Incorrect")) +  
theme(plot.title = element_text(size = 10, face = "bold", hjust = 0.5)) + 
theme(legend.position = "bottom", legend.text = element_text(size = 8, 
margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), 
                                                                                                                                 "lines")) + theme(legend.key.size = unit(0.4, "cm")) + ggtitle("(c) Correct and incorrect forecasts");
correct.map.e.2020

layout_matrix <- matrix(c(1, 1, 2, 2, 3, 3, 4, 4, 6, 5, 5, 6), nrow = 3, byrow = TRUE)
map.2020 <- grid.arrange(ev.p.g, ev.a.g, forecast.map.e.2020, outcome.map.e.2020, correct.map.e.2020, 
layout_matrix = layout_matrix); map.2020

map.2020

dev.off()

################################################################################
## 2024 ELECTION
################################################################################

# Merge with 2024 Biden forecast

forecasts_2024_Biden <- subset(forecasts_2024_Biden, select = c(fstatewinner_2024))

forecasts_2024_Biden <- forecasts_2024_Biden %>% 
  rename(fstatewinner_2024_Biden = fstatewinner_2024)

forecasts_2024 <- cbind(forecasts_2024_e,forecasts_2024_Biden)

# Figure 5. Two-Party Vote Share and Electoral College Vote Forecasts by State (Extended Model), 2024 (see article)

tiff("2024_map.tiff", units="in", width=7.5, height=7.5, res=1200)

harris.2024 <- statebins(forecasts_2024, value_col = "harris_v_2024", round = TRUE, name = "                   
                        Vote share  ", palette = "Blues", 
direction = 1) + theme_statebins() + theme(legend.position = "bottom") + 
theme(plot.margin = unit(c(0,0,0,0), "lines")) + ggtitle("(a) Kamala Harris predicted vote share") + 
theme(plot.title = element_text(size = 11, face = "bold"), plot.subtitle = element_text(size = 11)); harris.2024

trump.2024 <- statebins(forecasts_2024, value_col = "trump_v_2024", round = TRUE, name = "                   
                        Vote share  ", palette = "Reds", 
direction = 1) + theme_statebins() + theme(legend.position = "bottom") + 
theme(plot.margin = unit(c(0,0,0,0), "lines")) + ggtitle("(b) Donald Trump predicted vote share") +
theme(plot.title = element_text(size = 11, face = "bold"), plot.subtitle = element_text(size = 11)); trump.2024

forecast.winner.2024 <- mutate(forecasts_2024, value = fstatewinner_2024) %>% 
statebins(font_size = 3, dark_label = "white", light_label = "white", round = TRUE,
ggplot2_scale_function = scale_fill_manual, name = "                  ",
values = c("Kamala Harris" = "#2166ac", "Donald Trump" = "#b2182b")) + theme_statebins() + 
theme(legend.position = "bottom", legend.text = element_text(size = 11, 
margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), "lines")) + 
ggtitle("(c) Predicted winner by state: Harris v. Trump") + theme(plot.title = element_text(size = 11, 
face = "bold"), plot.subtitle = element_text(size = 11)) + theme(legend.key.size = unit(0.4, "cm")) + 
labs(subtitle = "     Trump: 341, Harris: 197"); forecast.winner.2024

forecast.winner.biden.2024 <- mutate(forecasts_2024, value = fstatewinner_2024_Biden) %>% statebins(font_size = 3,
dark_label = "white", light_label = "white", round = TRUE, ggplot2_scale_function = scale_fill_manual, name = "             ", 
values = c("Joe Biden" = "#2166ac", "Donald Trump" = "#b2182b")) + theme_statebins() + 
theme(legend.position = "bottom", legend.text = element_text(size = 11, 
margin = margin(r = 10, unit = "pt"))) + theme(plot.margin = unit(c(0,0,0,0), "lines")) + 
ggtitle("(d) Predicted winner by state: Biden v. Trump") + theme(plot.title = element_text(size = 11,
face = "bold"), plot.subtitle = element_text(size = 11)) + theme(legend.key.size = unit(0.4, "cm")) +
labs(subtitle = "     Trump: 216, Biden: 322"); forecast.winner.biden.2024

map.2024 <- ggarrange(harris.2024, trump.2024, NULL, NULL, forecast.winner.2024, forecast.winner.biden.2024, 
ncol = 2, nrow = 3, heights = c(1, 0.2, 1, 1, 0.2, 1)); map.2024

#dev.copy(pdf, "map_e_2024.pdf")
dev.off()

map.2024 <- ggarrange(ev.p.g, NULL, harris.2024, trump.2024, forecast.winner.2024, forecast.winner.biden.2024, 
                      ncol = 2, nrow = 3, heights = c(1, 1, 1, 0.2, 0.2, 1)); map.2024


################################################################################
## 2024 ELECTION LIKELIHOOD (BIDEN)
################################################################################

chance <- data.frame(
  gr = c("G", "G", "G", "G", "G", "G", "G"),
  party = c("1", "2", "3", "4", "5", "6", "7"),
  ev = c(166, 11, 68, 70, 33, 86, 104)
)

tiff("chance.tiff", units="in", width=8, height=7, res=1200)

chance.g <- ggplot(chance, aes(x = gr, y = ev)) + ggtitle("(b) Joe Biden v. Donald Trump") + 
  geom_col(aes(fill = party), show.legend=FALSE, width = .1) + 
  scale_fill_manual(values=c('#7F0000','#CC3333','#E88471','#FFEDA0',
                             '#D1E5F0', '#2B8CBE', '#084081'))

chance.g <- chance.g + coord_flip() + theme_void() + 
  theme(plot.margin = unit(c(0,0,-1.6,0), "inches"))

chance.g <- chance.g + annotation_custom(cursor, xmin=1, xmax=1.1, ymin=263, ymax=277)
chance.g <- chance.g + geom_text(aes(x = 1.08, y = 270, label = "270"), color = "black", size = 4.5)
chance.g <- chance.g + geom_text(aes(x = 1, y = 52, label = "104"), color = "white", size = 4)
chance.g <- chance.g + geom_text(aes(x = 1, y = 147, label = "86"), color = "white", size = 4)
chance.g <- chance.g + geom_text(aes(x = 1, y = 206, label = "33"), color = "black", size = 4)
chance.g <- chance.g + geom_text(aes(x = 1, y = 257, label = "70"), color = "black", size = 4)
chance.g <- chance.g + geom_text(aes(x = 1, y = 330, label = "68"), color = "black", size = 4)
chance.g <- chance.g + geom_text(aes(x = 1, y = 366, label = "11"), color = "white", size = 4)
chance.g <- chance.g + geom_text(aes(x = 1, y = 446, label = "166"), color = "white", size = 4)

chance.g

dev.off()


################################################################################
## 2024 ELECTION LIKELIHOOD (HARRIS)
################################################################################

chance.harris <- data.frame(
  gr = c("G", "G", "G", "G", "G", "G", "G"),
  party = c("1", "2", "3", "4", "5", "6", "7"),
  ev = c(252, 60, 24, 22, 0, 39, 141)
)

tiff("chance_harris.tiff", units="in", width=8, height=7, res=1200)

chance.harris.g <- ggplot(chance.harris, aes(x = gr, y = ev)) + ggtitle("Kamala Harris v. Donald Trump") + 
  geom_col(aes(fill = party), show.legend=FALSE, width = .1) + 
  scale_fill_manual(values=c('#7F0000','#CC3333','#E88471','#FFEDA0',
                             '#D1E5F0', '#2B8CBE', '#084081'))

chance.harris.g <- chance.harris.g + coord_flip() + theme_void() + 
  theme(plot.margin = unit(c(0,0,-1.6,0), "inches"))

chance.harris.g <- chance.harris.g + annotation_custom(cursor, xmin=1, xmax=1.1, ymin=263, ymax=277)
chance.harris.g <- chance.harris.g + geom_text(aes(x = 1.08, y = 270, label = "270"), color = "black", size = 4.5)
chance.harris.g <- chance.harris.g + geom_text(aes(x = 1, y = 70, label = "141"), color = "white", size = 4)
chance.harris.g <- chance.harris.g + geom_text(aes(x = 1, y = 160, label = "39"), color = "white", size = 4)
#chance.harris.g <- chance.harris.g + geom_text(aes(x = 1, y = 118.5, label = "35"), color = "black", size = 4)
chance.harris.g <- chance.harris.g + geom_text(aes(x = 1, y = 191, label = "22"), color = "black", size = 4)
chance.harris.g <- chance.harris.g + geom_text(aes(x = 1, y = 215, label = "24"), color = "black", size = 4)
chance.harris.g <- chance.harris.g + geom_text(aes(x = 1, y = 256, label = "60"), color = "white", size = 4)
chance.harris.g <- chance.harris.g + geom_text(aes(x = 1, y = 412, label = "252"), color = "white", size = 4)

chance.harris.g

dev.off()


################################################################################
## BEFORE-THE-FACT FORECASTS, 2004-2020: ALTERNATIVE MODEL
################################################################################

# Electoral College

df.ev2 <- data.frame(ev = c(373, 236, 200, 332, 206, 212, 
                            266, 286, 173, 332, 227, 232),
                     forecast = c("F", "F", "F", "F", "F", "F",
                                  "A", "A", "A", "A", "A", "A"),
                     year = c("2000", "2004", "2008", "2012", "2016", "2020")
)

tiff("ev_a_g.tiff", units="in", width=10, height=8.5, res=1200)

ev2.g <- df.ev2 %>% ggplot(aes(x=year, y=ev, fill=year)) + 
  geom_bar(stat="identity", alpha=.25, data = . %>% filter(forecast == "A"), width = 0.9, show.legend = FALSE) +
  geom_bar(stat="identity", data = . %>% filter(forecast == "F"), width = 0.6, show.legend = FALSE) +
  scale_y_continuous("Electoral Votes", limits=c(0,500), 
                     breaks=seq(0,500,100)) +
  ggtitle("   Multilevel Model") + geom_hline(yintercept=270, linetype="dashed") +
  scale_fill_manual(values=c("#0044C9", "#E81B23","#E81B23","#0044C9","#0044C9", "#E81B23")) + theme_minimal()

ev2.g <- ev2.g + 
  theme(axis.title.y=element_blank(),
        axis.title.x=element_text(margin=margin(t=8, r=0, b=0, l=0), size = 14),
        axis.text.x=element_text(size = 14),
        axis.text.y=element_text(size = 14, color="black"),
        plot.title = element_text(size = 15, margin=margin(t=0, r=0, b=8, l=0))) + coord_flip()

ev2.g <- ev2.g + annotation_custom(gore, xmin=.6, xmax=1.4, ymin=6, ymax=42)
ev2.g <- ev2.g + annotation_custom(bush, xmin=1.6, xmax=2.4, ymin=6, ymax=42)
ev2.g <- ev2.g + annotation_custom(mccain, xmin=2.6, xmax=3.4, ymin=6, ymax=42)
ev2.g <- ev2.g + annotation_custom(obama, xmin=3.6, xmax=4.4, ymin=6, ymax=42)
ev2.g <- ev2.g + annotation_custom(clinton, xmin=4.6, xmax=5.4, ymin=6, ymax=42)
ev2.g <- ev2.g + annotation_custom(trump2, xmin=5.6, xmax=6.4, ymin=6, ymax=42)
ev2.g <- ev2.g + geom_text(aes(x = 6.45, y = 305, label = "270 to win"), color = "black", size = 4.5)
ev2.g <- ev2.g + geom_text(aes(x = 6.21, y = 50, label = "212"), color = "white", size = 4)
ev2.g <- ev2.g + geom_text(aes(x = 6.37, y = 50, label = "232"), color = "black", size = 4)
ev2.g <- ev2.g + geom_text(aes(x = 5.21, y = 50, label = "206"), color = "white", size = 4)
ev2.g <- ev2.g + geom_text(aes(x = 5.37, y = 50, label = "227"), color = "black", size = 4)
ev2.g <- ev2.g + geom_text(aes(x = 4.21, y = 50, label = "332"), color = "white", size = 4)
ev2.g <- ev2.g + geom_text(aes(x = 4.37, y = 50, label = "332"), color = "black", size = 4)
ev2.g <- ev2.g + geom_text(aes(x = 3.21, y = 50, label = "200"), color = "white", size = 4)
ev2.g <- ev2.g + geom_text(aes(x = 3.37, y = 50, label = "173"), color = "black", size = 4)
ev2.g <- ev2.g + geom_text(aes(x = 2.21, y = 50, label = "236"), color = "white", size = 4)
ev2.g <- ev2.g + geom_text(aes(x = 2.37, y = 50, label = "286"), color = "black", size = 4)
ev2.g <- ev2.g + geom_text(aes(x = 1.21, y = 50, label = "373"), color = "white", size = 4)
ev2.g <- ev2.g + geom_text(aes(x = 1.37, y = 50, label = "266"), color = "black", size = 4)

ev2.g

dev.off()